Mega-analysis of the Interoceptive Accuracy Scale (IAS)

Study 1

library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)

Data Preparation

Sample 1

• Murphy (2020)
  • https://osf.io/3m5nh

df1a <- haven::read_sav("../data/Murphy2020/Study 1.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

df1b <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

Sample 2

• Gaggero (2021)
  • https://osf.io/5x9sg

load("../data/Gaggero2021/data.rda")
df2 <- data |> 
  mutate(Gender = as.character(gender)) |> 
  mutate(across(starts_with("IAS "), as.numeric)) |>
  rename(
    Age=age, Heart = `IAS 1`, Hungry = `IAS 2`, Breathing = `IAS 3`, Thirsty = `IAS 4`,
    Urinate = `IAS 5`, Defecate = `IAS 6`, Taste = `IAS 7`, Vomit = `IAS 8`, Sneeze = `IAS 9`,
    Cough = `IAS 10`, Temp = `IAS 11`, Sex_arousal = `IAS 12`, Wind = `IAS 13`, Burp = `IAS 14`,
    Muscles = `IAS 15`, Bruise = `IAS 16`, Pain = `IAS 17`, Blood_Sugar = `IAS 18`,
    Affective_touch = `IAS 19`, Tickle = `IAS 20`, Itch = `IAS 21`)
rm(data)

Sample 3

• Campos (2022) - Study 1
  • https://osf.io/j6ef3

df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )

Sample 4

• Todd (2022)
  • https://osf.io/ms354/

df4 <- haven::read_sav("../data/Todd2022/CompleteData_new.sav") |>
  rename_with(\(x) str_remove(x, "IAS_"), .cols = starts_with("IAS_")) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )

Sample 5

• Arslanova (2022)
  • https://osf.io/mp3cy/
  • note: Final.xlsx includes information about which participants passed the attention checks

df5_attention <- readxl::read_excel("../data/Arslanova2022/Participants2.xlsx", 
    sheet = "FINAL") |>
  filter(ActorHR == 1) |>
  select(Subj.ID, Gender, Age) |>
  mutate(Gender = as.numeric(Gender), Age = as.numeric(Age)) |>
  filter(!is.na(Age)) |> #error on the documentation of this participant - attention check failed but not noted
  mutate(Gender = case_when(
    Gender== 0 ~ "Male",
    Gender== 1 ~ "Female" # based on paper reporting 65 women
  ))

New names:
• `Why No?` -> `Why No?...4`
• `Why No?` -> `Why No?...6`
• `` -> `...12`
• `` -> `...18`
• `` -> `...19`
• `` -> `...20`
• `` -> `...21`
• `` -> `...22`
• `` -> `...23`
• `` -> `...24`
• `` -> `...25`
• `` -> `...26`
• `` -> `...27`
• `` -> `...28`

df5 <- readxl::read_excel("../data/Arslanova2022/Murphy_Iacc_scale.xlsx")  |>
  filter(str_detect(Question.Key, "IAC_")) |> 
  filter(str_detect(Question.Key, "-quantised")) |> 
  pivot_wider(names_from = Question.Key, values_from = Response) |>
  rename_with(\(x) str_remove(x, "IAC_"), .cols = starts_with("IAC_")) |> 
  rename_with(\(x) str_remove(x, "-quantised"), .cols = everything()) |> 
  rename(Sex_arousal = SexuallyAroused, Itch = Itchy, Sex_arousal = SexuallyAroused,
         Temp = HotCold, Tickle = Ticklish, Breathing= BreatheFast, 
         Affective_touch= PleasantAffectionate, Blood_Sugar= BloodSugar, 
         Muscles=TiredMuscles, Heart= FastHeart, Taste=Tastes) |> 
  rename(Subj.ID = "Participant.Public.ID") |>
  select(-starts_with("Participant")) 

df5 <- df5 |>
  filter(Subj.ID %in% df5_attention$Subj.ID) |>
  select(-Subj.ID) |>
  mutate(across(everything(), as.numeric))

Sample 6

• Brand (2022)
  • https://osf.io/xwz6g/

df6 <- haven::read_sav("../data/Brand2022/LatentVariableApproach.sav") |> 
  select(-SERIAL) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

Sample 7

• Brand (2023)
  • https://osf.io/3f2h6/

df7a <- haven::read_sav("../data/Brand2023/Data_Giessen.sav") |> 
  select(-COUNTRY_OTHER) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(SEX == 1, "Male", ifelse(SEX == 2, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

df7b <- haven::read_sav("../data/Brand2023/Data_Mainz.sav") |> 
  select(-COUNTRY_OTHER, -EDUCATION_OTHER, -Sample) |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )
  
  
df7c <- haven::read_sav("../data/Brand2023/Data_PotVie.sav") |> 
  select(-ID) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IA02_01, Hungry = IA02_02, Breathing = IA02_03, Thirsty = IA02_04,
    Urinate = IA02_05, Defecate = IA02_06, Taste = IA02_07, Vomit = IA02_08, Sneeze = IA02_09,
    Cough = IA02_10, Temp = IA02_11, Sex_arousal = IA02_12, Wind = IA02_13, Burp = IA02_14,
    Muscles = IA02_15, Bruise = IA02_16, Pain = IA02_17, Blood_Sugar = IA02_18,
    Affective_touch = IA02_19, Tickle = IA02_20, Itch = IA02_21
  )

Sample 8

• Lin (2023)
  • https://osf.io/3eztd
  • Note: No tickle because same Chinese character

df8a <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |>
  select(-sex) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |>
  rename(
    Age = age,
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH
  ) 


df8b <- haven::read_sav("../data/Lin2023/Study 2.sav") |>
  mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |>
  rename(
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH)

Sample 9

• VonMohr (2023) - Study 3
  • https://osf.io/7p9u5/

df9 <- read.csv("../data/VonMohr2023/DataSet_study3.csv")
df9 <- df9[complete.cases(select(df9, starts_with("IAS_"))), ]
df9 <- df9 |>
  mutate(Gender = as.character(ifelse(GENDER == "Man", "Male", ifelse(GENDER == "Woman", "Female", "Other")))) |>
  rename(
    Age=AGE, Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

Sample 10

• Makowski
  • Note: First sample has some missing items (No Taste, Cough, Blood_Sugar)
  • Note: Analog scales
  • note df10c = primals data

df10a <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  rename(
    Gender = Sex, Heart = Item_IAS_Interoception_1, Hungry = Item_IAS_Interoception_2,
    Breathing = Item_IAS_Interoception_3, Thirsty = Item_IAS_Interoception_4,
    Temp = Item_IAS_Interoception_5, Sex_arousal = Item_IAS_Interoception_6,
    Urinate = Item_IAS_Elimination_1, Defecate = Item_IAS_Elimination_2,
    Vomit = Item_IAS_Elimination_3, Wind = Item_IAS_Expulsion_1,
    Burp = Item_IAS_Expulsion_2, Sneeze = Item_IAS_Expulsion_3,
    Muscles = Item_IAS_Nociception_1, Bruise = Item_IAS_Nociception_2,
    Pain = Item_IAS_Nociception_3, Affective_touch = Item_IAS_Skin_1,
    Tickle = Item_IAS_Skin_2, Itch = Item_IAS_Skin_3
  ) |>
  filter(!is.na(Urinate))

df10b <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
  rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21,
    typeSample = Sample
  )

df10c <- read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionPrimals/refs/heads/main/data/data_participants.csv") |>
  select("Participant" = "participant_id", "Gender", "Age", "Ethnicity", starts_with(c("IAS","MAIA","PI"))) |>
   rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  ) |>
  filter(!is.na(Heart)) # participant is outlier and had total IAS scores below 0.4 

df10d <- read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionScale/refs/heads/main/study2/data/rawdata_participants.csv") |>
  select("Participant", "Gender", "Age", "Ethnicity", starts_with(c("IAS","MAIA","PI", "PHQ4"))) |>
   rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  ) |> 
  filter(Participant %in% read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionScale/refs/heads/main/study2/data/data_participants.csv")$Participant)

# pi_vars <- c("PI_Enticing", "PI_Alive", "PI_Safe", "PI_Good", "PI_Changing", 
#              "PI_Hierarchical", "PI_Understandable", "PI_AboutMe", "PI_Abundant", 
#              "PI_Acceptable", "PI_Beautiful", "PI_Cooperative", "PI_Funny", 
#              "PI_Harmless", "PI_Improvable", "PI_Intentional", "PI_Interconnected", 
#              "PI_Interesting", "PI_Just", "PI_Meaningful", "PI_NeedsMe", 
#              "PI_Pleasurable", "PI_Progressing", "PI_Regenerative", "PI_Stable", 
#              "PI_WorthExploring")
# 
# 
# correlation::correlation(
#   df10c[, "IAS_Total", drop = FALSE],  
#   data2 = select(df10c, all_of(pi_vars)), p_adjust = "none"
# ) |> 
#   # Formatting correlation results
#   mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
#   
#   # Releveling factors for visualization
#   mutate(Parameter2 = fct_relevel(Parameter2, rev(pi_vars)),
#          Parameter1 = fct_relevel(Parameter1, "IAS_Total")) |>
#   
#   # Plotting the correlation matrix
#   ggplot(aes(x = Parameter1, y = Parameter2)) +
#   geom_tile(aes(fill = r), color = "white") +  
#   geom_text(aes(label = val), size = 3) +  
#   labs(title = "Correlation Matrix") +
#   scale_fill_gradient2(low = "blue", high = "red", mid = "white", limit = c(-1, 1), midpoint = 0, guide = guide_colourbar(ticks = FALSE)) +  
#   theme_minimal() +
#   theme(
#     legend.title = element_blank(),
#     axis.title.x = element_blank(),
#     axis.title.y = element_blank(),
#     axis.text.x = element_text(hjust = 1)
#   )

Sample 11

• Poreiro et al.,
  • Study 1: https://osf.io/49wbv
  • Study 2: unpublished

df11a <- haven::read_sav("../data/Giulia/Interoceptive Attention ESM.sav") |> 
  select(AGE, GENDER, contains("IAS_ACC"), -IAS_ACC_20, starts_with(c("MAIA", "TAS20", "DEP", "ANX")),  -ANX_21, -TAS20_21) |> #anx_21 , tas20_21 are attention checks
  mutate_all(as.numeric) |>
  rename(
    Age = AGE, Gender = GENDER,
    Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4,
    Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7, Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9,
    Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
    Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18,
    Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_21, Itch = IAS_ACC_22
  ) |>
   mutate(Gender = case_when(
    Gender == 1 ~ "Male",
    Gender == 2 ~ "Female",
    Gender %in% c(3, 5) ~ "Other",
    Gender == 4 ~ "NA"
  )) |> ## compute dimensions for the TAS
  mutate(TAS_DIF = TAS20_1 + TAS20_3 + TAS20_6 + TAS20_7 + TAS20_9 + TAS20_13 + TAS20_14,
         TAS_DDF = TAS20_2 + (1-TAS20_4) + TAS20_11 + TAS20_12 + TAS20_17,
         TAS_EOT = (1-TAS20_5) + TAS20_8 + (1-TAS20_10) + TAS20_15 + TAS20_16 + (1-TAS20_18) + (1-TAS20_19) + TAS20_20)


df11b <- haven::read_sav("../data/Giulia/Sleep and Exteroceptive Interoceptive Sensitivity.sav") |> # 165 participants
  filter(Total_Attention_Fail == 0) |> # 32 participants removed based on failing one or more checks
  select(AGE, GENDER, contains("IAS_ACC"), ) |> 
  mutate_all(as.numeric) |> 
  rename( Age = AGE, Gender = GENDER, 
          Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4, Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7,
          Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9, Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
          Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18, Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_20, 
          Itch = IAS_ACC_21 )  |> 
  mutate(Gender = case_when(
    Gender == 1 ~ "Male", 
    Gender == 2 ~"Female", 
    Gender %in% c(3, 5) ~ "Other", 
    Gender == 4 ~ "NA", )) |>
  na.omit()


# Save all
data <- list(df1a=df1a, df1b=df1b, df2=df2, df3=df3, df4=df4, df5=df5, df6=df6, df7a=df7a, df7b=df7b, df7c=df7c, df8a=df8a, df8b=df8b, df9=df9, df10a=df10a, df10b=df10b, df10c=df10c, df11a=df11a, df11b=df11b)
save(data, file = "../data/data.RData")

Participants

• Sample 1a: Data from Murphy’s (2020) study 1, downloaded from OSF, included 451 participants (Mean age = 25.8, SD = 8.4, range: [18, 69]; Gender: 69.4% women, 29.5% men, 1.11% non-binary).
• Sample 1b: Data from Murphy’s (2020) study 6, downloaded from OSF, included 375 participants (Mean age = 35.3, SD = 16.9, range: [18, 91]; Gender: 70.1% women, 28.5% men, 1.33% non-binary).
• Sample 2: Data from Gaggero’s(2020) study, downloaded from OSF, included 814 participants (Mean age = 24.9, SD = 5.3, range: [18, 58], 0.2% missing; Gender: 60.3% women, 39.4% men, 0.25% non-binary).
• Sample 3: Data from Campos’s(2022) study, downloaded from OSF, included 515 participants (Mean age = 30.7, SD = 10.5, range: [18, 72]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 59.6% women, 40.4% men, 0.00% non-binary).
• Sample 4: Data from Todd’s(2022) study, downloaded from OSF, included 802 participants ().
• Sample 5: Data from Arslanova (2022) study, downloaded from OSF, included 143 participants (Mean age = 28.5, SD = 7.6, range: [18, 73]; Gender: 45.5% women, 54.5% men, 0.00% non-binary).
• Sample 6: Data from Brand’s(2022) study, downloaded from OSF, included 619 participants (Mean age = 43.9, SD = 14.5, range: [18, 78]; Gender: 78.7% women, 20.2% men, 1.13% non-binary).
• Sample 7a: Data from Brand’s(2023) study, downloaded from OSF, included 522 participants (Mean age = 23.4, SD = 6.7, range: [18, 79]; Gender: 79.5% women, 19.7% men, 0.77% non-binary).
• Sample 7b: Data from Brand’s(2023) study, downloaded from OSF, included 1993 participants (Mean age = 32.0, SD = 12.6, range: [16, 81]; Gender: 77.7% women, 21.7% men, 0.60% non-binary).
• Sample 7c: Data from Brand’s(2023) study, downloaded from OSF, included 808 participants (Mean age = 27.3, SD = 9.3, range: [18, 72], 0.5% missing; Gender: 68.9% women, 30.2% men, 0.87% non-binary).
• Sample 8a: Data from Lin’s(2023) study, downloaded from OSF, included 1166 participants (Mean age = 32.5, SD = 8.4, range: [16, 60]; Gender: 57.0% women, 43.0% men, 0.00% non-binary).
• Sample 8b: Data from Lin’s(2023) study, downloaded from OSF, included 500 participants (Mean age = 37.4, SD = 7.4, range: [20, 60]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 56.2% women, 43.8% men, 0.00% non-binary).
• Sample 9: Data from VonMohr’s(2023) study 3, downloaded from OSF, included 21843 participants (Mean age = 56.5, SD = 14.4, range: [18, 93], 0.2% missing; Gender: 73.2% women, 25.1% men, 1.55% non-binary, 0.15% missing).
• Sample 10a: Data from [Makowski’s(2023)] study , downloaded from OSF, included 485 participants (Mean age = 30.1, SD = 10.1, range: [18, 73]; Gender: 50.3% women, 49.7% men, 0.00% non-binary; Education: Bachelor, 45.15%; Doctorate, 1.86%; High school, 34.43%; Master, 15.88%; Other, 2.47%; Prefer not to say, 0.21%).
• Sample 10b: Data from [Makowski’s(2023)] study , downloaded from OSF, included 836 participants (Mean age = 25.1, SD = 11.3, range: [18, 76], 0.1% missing; Gender: 73.8% women, 22.6% men, 2.87% non-binary, 0.72% missing; Education: Bachelor, 22.85%; Doctorate, 2.15%; High School, 63.52%; Master, 6.22%; missing, 0.24%; Other, 5.02%).
• Sample 10c: Data from [Makowski’s(2023)] study , downloaded from OSF, included 146 participants (Mean age = 21.1, SD = 4.3, range: [18, 50], 0.7% missing; Gender: 80.8% women, 15.1% men, 2.74% non-binary, 1.37% missing).
• Sample 11a: Data from Poreiro’s(2024) study , included 107 participants (Mean age = 26.8, SD = 9.2, range: [18, 57]; Gender: 74.8% women, 23.4% men, 1.87% non-binary)
• Sample 11b: Data from [Poreiro’s] study , included 131 participants (Mean age = 30.9, SD = 12.0, range: [18, 60]; Gender: 76.3% women, 22.9% men, 0.76% non-binary)

Total N = 32256.

desc <- function(df){
  summary <- df |> summarise(
      n = n(),
      mean = mean(Age, na.rm = TRUE),
      sd = sd(Age, na.rm = TRUE),
      min = min(Age, na.rm = TRUE),
      max = max(Age, na.rm = TRUE),
      female = sum(Gender == "Female", na.rm = TRUE)
    )
    return(summary)
}

desc_df1a <- desc(df1a)
desc_df1b <- desc(df1b)
desc_df2 <- desc(df2)
desc_df3 <- desc(df3)
#desc_df4 <- desc(df4) # no demographic data available
desc_df5 <- desc(df5_attention)
desc_df6 <- desc(df6)
desc_df7a <- desc(df7a)
desc_df7b <- desc(df7b)
desc_df7c <- desc(df7c)
desc_df8a <- desc(df8a)
desc_df8b <- desc(df8b)
desc_df9 <- desc(df9)
desc_df10a <- desc(df10a)
desc_df10b <- desc(df10b)
desc_df10c <- desc(df10c)
desc_df11a <- desc(df11a)
desc_df11b <- desc(df11b)

desc_total <- rbind(desc_df1a, desc_df1b, desc_df2, desc_df3, desc_df5, desc_df6, desc_df7a, desc_df7b, desc_df7c, desc_df8a, desc_df8b, desc_df9, desc_df10a, desc_df10b, desc_df10c, desc_df11a, desc_df11b)

## calculate weighted mean

total_mean = sum(desc_total$mean * desc_total$n) / sum(desc_total$n)
total_sd = sum(desc_total$sd * desc_total$n) / sum(desc_total$n)
perce_female = sum(desc_total$female)/sum(desc_total$n) *100

library(gt)

# APA style ####
gt_apastyle <- function(gt_table, font.size=12) {
  gt_table  |> 
    gt::opt_table_lines(extent = "none") %>%
    gt::tab_options(
      heading.border.bottom.width = 2,
      heading.border.bottom.color = "black",
      heading.border.bottom.style = "solid",
      table.border.top.color = "black",
      table.border.top.style = "solid",
      table.border.top.width = 2,  
      table_body.hlines.color = "white",
      table_body.border.top.color = "black",
      table_body.border.top.style = "solid",
      table_body.border.top.width = 2,
      heading.title.font.size = font.size,
      table.font.size = font.size,
      heading.subtitle.font.size = font.size,
      table_body.border.bottom.color = "black",
      table_body.border.bottom.width = 2,
      table_body.border.bottom.style = "solid",
      column_labels.border.bottom.color = "black",
      column_labels.border.bottom.style = "solid",
      column_labels.border.bottom.width = 1,
      latex.use_longtable = FALSE
    ) |> 
      gt::opt_table_font(font = "times")
}

make_age <- function(age) {
  age <- as.numeric(age) 
  
  mean(age, na.rm = TRUE) |> 
    insight::format_value(digits=1) |> 
    paste0(" ± ", 
           insight::format_value(sd(age, na.rm = TRUE), digits=1)) 
} 

# Create the dataframe
table <- data.frame(
  Reference = c('Murphy et al., (2020)', '', '', 'Gaggero et al., (2021)', 'Campos et al., (2022)', 'Todd et al., (2022)', 'Arslanova et al., (2022)', 'Brand et al., (2022)', 'Brand et al., (2023)', '', '', '', 'Lin et al., (2023)', '', '', 'VonMohr et al., (2023)', 'Makowski et al., (2023a)', 'Makowski et al., (2023b)', 'Makowski et al., (2023c)', 'Poerio et al., (2024)', 'Poerio et al., unpublished', 'Total'),
  Sample = c('', 'Sample 1a', 'Sample 1b', 'Sample 2', 'Sample 3', 'Sample 4', 'Sample 5', 'Sample 6', '', 'Sample 7a', 'Sample 7b', 'Sample 7c', '', 'Sample 8a', 'Sample 8b', 'Sample 9', 'Sample 10', 'Sample 11', 'Sample 12', 'Sample 13', 'Sample 14', 'Sample 15'),   
  Language = c('', 'English', 'English' , 'English and Italian', 'Portuguese', 'English', 'English', 'German', '', 'German', 'German', 'German', '', 'Chinese', 'Chinese', 'English', 'English', 'English', 'English', 'English', 'English',''),
  N = c('', nrow(df1a), nrow(df1b), nrow(df2), nrow(df3), nrow(df4), nrow(df5_attention), nrow(df6),'', nrow(df7a), nrow(df7b), 802,'', nrow(df8a), nrow(df8b), nrow(df9), nrow(df10a), nrow(df10b), nrow(df10c), nrow(df11a), nrow(df11b), 32214),
  'Difference' = c('','', '', '', '', '', '', '', '', '', '', '', '', 'Collapsed "Itch" and "Tickling"', 'Collapsed "Itch" and "Tickling"', '', 'Analog scales', 'Analog scales', 'Analog scales', '', '', ''),
  Age = c('', make_age(df1a$Age), make_age(df1b$Age), make_age(df2$Age), make_age(df3$Age), "48.6.6 ± 14.1*", make_age(df5_attention$Age), make_age(df6$Age),'', make_age(df7a$Age), make_age(df7b$Age), make_age(df7c$Age),'', make_age(df8a$Age), make_age(df8b$Age), make_age(df9$Age), make_age(df10a$Age), make_age(df10b$Age), make_age(df10c$Age), make_age(df11a$Age), make_age(df11b$Age), '48.6 ± 13.1'),
  Range = c('', '18-69', '18-91', '18-58', '18-72', '18-92*', '18-73', '18-78', '', '18-79', '16-81', '18-72','', '16-60', '20-60', '18-93', '18-73', '17-76','18-50', '18-57', '18-60', '17-93'),
  Female_Percentage = c('', '69.4%', '70.1%', '60.3%', '59.6%', '50%*', '46.8%', '78.7%', '', '79.5%', '77.7%', '68.9%', '', '57.0%', '56.2%', '73.2%', '50.3%', '53.0%', '76%', '74.8%', '75.9%', '71.6%'),
  Availability= c('osf.io/3m5nh', '', '', 'osf.io/5x9sg', 'osf.io/j6ef3', 'osf.io/ms354', 'osf.io/mp3cy', 'osf.io/xwz6g', 'osf.io/3f2h6', '', '', '','osf.io/3eztd','', '', 'osf.io/7p9u5', 'github.com/RealityBending/IllusionGameReliability', 'github.com/DominiqueMakowski/PHQ4R', 'github.com/RealityBending/InteroceptionPrimals', 'osf.io/49wbv', '','')
) 

table_apa <- table |> 
  gt() |>
  cols_align(align = c("right"), columns = "Age") |> 
  cols_label(Age = "Age (Mean  ± SD)", Female_Percentage = "Female %") |> 
  tab_footnote("* Information taken from the sample description of relevant paper rather than recomputed.")
  
gt_apastyle(table_apa, font.size=12)
gtsave(gt_apastyle(table_apa, font.size=9), "figures/table1.tex")

# gtsave(table_apa, "figures/table1.png")

Item Selection

Distribution

The items with the differing distribution pattern (with a lower mode around 2/5) are “Affective Touch”, “Blood Sugar” and “Bruise”.

vars <- c(
  "Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
  "Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch"
)

dens1a <- estimate_density(select(df1a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1a")
dens1b <- estimate_density(select(df1b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1b")
dens2 <- estimate_density(select(df2, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 6")
dens7a <- estimate_density(select(df7a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7a")
dens7b <- estimate_density(select(df7b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7b")
dens7c <- estimate_density(select(df7c, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7c")
dens8a <- estimate_density(select(df8a, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8a")
dens8b <- estimate_density(select(df8b, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8b")
dens9 <- estimate_density(select(df9, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 9")
dens10a <- estimate_density(select(df10a, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method = "kernSmooth") |>
  mutate(Sample = "Sample 10a",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10b <- estimate_density(select(df10b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 10b",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10c <- estimate_density(select(df10c, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 10c",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens11a <- estimate_density(select(df11a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 11a")
dens11b <- estimate_density(select(df11b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 11b")


p1_data <- rbind(dens1a, dens1b, dens2, dens3, dens4, dens5, dens6, dens7a, dens7b, dens7c, dens8a, dens8b, dens9, dens10a, dens10b, dens10c, dens11a, dens11b) 

linetype <- setNames(rep("solid", length(vars)), vars)
linetype["Affective touch"] <- "dotted"
linetype["Blood sugar"] <- "dashed"
linetype["Sex arousal"] <- "solid"
linetype["Bruise"] <- "dashed"

# proportion of non-1 and non-5 values 
p1_data |>
  filter(x > 1 & x < 5) |>  # Keep only x values strictly between 1 and 5
  summarise(proportion = sum(y) / sum(p1_data$y))  # Normalize by total density

  proportion
1  0.9994812

p1 <- p1_data |>
  mutate(Sample = fct_relevel(Sample, "Sample 1a", "Sample 1b", "Sample 2", "Sample 3", "Sample 4", "Sample 5", "Sample 6", "Sample 7a", "Sample 7b", "Sample 7c", "Sample 8a", "Sample 8b", "Sample 9", "Sample 10a", "Sample 10b", "Sample 10c", "Sample 11a", "Sample 11b"),
         Parameter = ifelse(Parameter == "Affective_touch", "Affective touch", Parameter),
         Parameter = ifelse(Parameter == "Blood_Sugar", "Blood sugar", Parameter),
         Parameter = ifelse(Parameter == "Sex_arousal", "Sex arousal", Parameter)) |>
  # mutate(Dashed = ifelse(Parameter %in% c("Affective_touch", "Blood_Sugar", "Bruise"), TRUE, FALSE)) |> 
  ggplot(aes(x = x, y = y, color = Parameter)) +
  geom_line(aes(linetype = Parameter), linewidth = 0.7) +
  scale_color_material_d() +
  scale_linetype_manual(values = linetype)  +
  labs(x = "Response", y = "Distribution", color = "Item", linetype = "Item", title = "Item Distribution") +
  guides(color = guide_legend(ncol = 1)) +
  facet_wrap(~Sample, scales = "free_y", nrow = 3) +
  theme_minimal() +
  theme(axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        plot.title = element_text(face="bold"))

p1

data1a <- normalize(select(df1a, all_of(dens1a$Parameter)), range = c(1, 5))
data1b <- normalize(select(df1b, all_of(dens1b$Parameter)), range = c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range = c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range = c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range = c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range = c(1, 5))
data6 <- normalize(select(df6, all_of(dens6$Parameter)), range = c(1, 5))
data7a <- normalize(select(df7a, all_of(dens7a$Parameter)), range = c(1, 5))
data7b <- normalize(select(df7b, all_of(dens7b$Parameter)), range = c(1, 5))
data7c <- normalize(select(df7c, all_of(dens7c$Parameter)), range = c(1, 5))
data8a <- normalize(select(df8a, all_of(dens8a$Parameter)), range = c(1, 5))
data8b <- normalize(select(df8b, all_of(dens8b$Parameter)), range = c(1, 5))
data9 <- normalize(select(df9, all_of(dens9$Parameter)), range = c(1, 5))
data10a <- select(df10a, all_of(dens10a$Parameter))
data10b <- select(df10b, all_of(dens10b$Parameter))
data10c <- select(df10c, all_of(dens10c$Parameter))
data11a <- normalize(select(df11a, all_of(dens11a$Parameter)), range = c(1, 5))
data11b <- na.omit(normalize(select(df11b, all_of(dens11b$Parameter)), range = c(1, 5))) # remove NAS
row.names(data11b) <- NULL



data_all <- rbind(
  data1a, data1b, data2, data3, data4, data5, data6, data7a, data7b, data7c,
  mutate(data8a, Tickle = NA), 
  mutate(data8b, Tickle = NA), data9,
  mutate(data10a, Taste = NA, Cough = NA, Blood_Sugar = NA),  
  data10b,  data10c, data11a, data11b
) 

Correlations

An overall positive intercorrelation pattern, with no clear structure emerging.

All data

make_correlation <- function(df) {
  correlation::correlation(df, redundant = TRUE) |>
    correlation::cor_sort() |>
    # correlation::cor_lower() |>
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE)),
           r = ifelse(Parameter1 == Parameter2, NA, r),
           lbl = ifelse(Parameter1 == Parameter2, "", format_value(r, zap_small = TRUE, lead_zero = FALSE)),
           lbl = ifelse(p > .05, "", lbl),
           Param1 = fct_relevel(str_replace(Parameter1, "_", " "), str_replace(levels(Parameter1), "_", " ")),
           Param2 = fct_relevel(str_replace(Parameter2, "_", " "), str_replace(levels(Parameter2), "_", " "))) |>
    # mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x = Param1, y = Param2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = lbl), size = 3) +
    labs(title = "Correlation Matrix", subtitle = paste0("N = ", nrow(df))) +
    scale_fill_gradientn(
      colors = c("white", "#FFF9C4", "#FFF59D", "#FFEB3B", "#FFCA28", "#FF9800", "#FF5722", "#F44336", "#E91E63", "#9C27B0", "#673AB7", "#512DA8", "#311B92"),
      na.value = "#311B92",
      limits = c(0, 0.8),
      guide = guide_colorbar(ticks.colour = NA)
    ) +
    # scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks = FALSE)) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1),
      plot.title = element_text(face = "bold")
    )
}

p2 <- make_correlation(data_all)
p2

fig1 <- p1 / p2
ggsave("figures/Figure1.png", fig1, width=10, height=12, dpi=300, bg="white")

Sample 1a

make_correlation(data1a)

Sample 1b

make_correlation(data1b)

Sample 2

make_correlation(data2)

Sample 3

make_correlation(data3)

Sample 4

make_correlation(data4)

Sample 5

make_correlation(data5)

Sample 6

make_correlation(data6)

Sample 7a

make_correlation(data7a)

Sample 7b

make_correlation(data7b)

Sample 7c

make_correlation(data7c)

Sample 8a

make_correlation(data8a)

Sample 8b

make_correlation(data8b)

Sample 9

make_correlation(data9)

Sample 10a

make_correlation(data10a)

Sample 10b

make_correlation(data10b)

Sample 10c

make_correlation(data10c)

Sample 11a

make_correlation(data11a)

Sample 11b

make_correlation(data11b)

Unique Variable Analysis (UVA)

Unique Variable Analysis (Christensen, Garrido, & Golino, 2023) uses the weighted topological overlap measure (Nowick et al., 2009) on an estimated network. Values greater than 0.25 are determined to have considerable local dependence (i.e., redundancy) that should be handled (variables with the highest maximum weighted topological overlap to all other variables (other than the one it is redundant with) should be removed).

All data

uva0 <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva0

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.363

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.290
 Urinate Defecate 0.265

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Sneeze     Cough 0.230
  Heart Breathing 0.225
 Hungry   Thirsty 0.217

uva0$keep_remove

$keep
[1] "Itch"

$remove
[1] "Tickle"

Sample 1a

uva1a <- EGAnet::UVA(data = data1a, cut.off = 0.3)
uva1a

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.286
   Wind   Burp 0.275

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.244
 Urinate Defecate 0.241

uva1a$keep_remove

NULL

Sample 1b

uva1b <- EGAnet::UVA(data = data1b, cut.off = 0.3)
uva1b

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.350
 Sneeze  Cough 0.317
   Wind   Burp 0.309

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva1b$keep_remove

$keep
[1] "Cough" "Burp"  "Itch" 

$remove
[1] "Sneeze" "Wind"   "Tickle"

Sample 2

uva2 <- EGAnet::UVA(data = data2, cut.off = 0.3)
uva2

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
 Tickle      Itch 0.283
  Heart Breathing 0.252

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.230
 Urinate Defecate 0.223
  Hungry  Thirsty 0.200

uva2$keep_remove

NULL

Sample 3

uva3 <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva3

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.445
 Sneeze  Cough 0.318

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
   Wind   Burp 0.256

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Bruise Blood_Sugar 0.219
 Urinate    Defecate 0.217
   Heart   Breathing 0.208

uva3$keep_remove

$keep
[1] "Sneeze" "Itch"  

$remove
[1] "Cough"  "Tickle"

Sample 4

uva4 <- EGAnet::UVA(data = data4, cut.off = 0.3)
uva4

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.415
 Tickle      Itch 0.351

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
  Hungry  Thirsty 0.293
 Urinate Defecate 0.282
    Wind     Burp 0.275
  Sneeze    Cough 0.251

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i node_j   wto
 Muscles   Pain 0.205

uva4$keep_remove

$keep
[1] "Heart" "Itch" 

$remove
[1] "Breathing" "Tickle"   

Sample 5

uva5 <- EGAnet::UVA(data = data5, cut.off = 0.3)

Warning: Some variables did not belong to a dimension: Blood_Sugar, Affective_touch 

Use caution: These variables have been removed from the TEFI calculation

uva5

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

    node_i node_j   wto
    Tickle   Itch 0.247
    Bruise   Itch 0.246
 Breathing  Vomit 0.207

uva5$keep_remove

NULL

Sample 6

uva6 <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva6

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.416
  Sneeze    Cough 0.332
    Wind     Burp 0.314
  Tickle     Itch 0.303

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i  node_j   wto
 Hungry Thirsty 0.274

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.237

uva6$keep_remove

$keep
[1] "Urinate" "Sneeze"  "Wind"    "Itch"   

$remove
[1] "Defecate" "Cough"    "Burp"     "Tickle"  

Sample 7a

uva7a <- EGAnet::UVA(data = data7a, cut.off = 0.3)
uva7a

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Sneeze  Cough 0.434
 Tickle   Itch 0.376

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.295
 Urinate Defecate 0.294

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva7a$keep_remove

$keep
[1] "Sneeze" "Tickle"

$remove
[1] "Cough" "Itch" 

Sample 7b

uva7b <- EGAnet::UVA(data = data7b, cut.off = 0.3)
uva7b

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.317
  Tickle     Itch 0.311

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.281
 Sneeze     Cough 0.261

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.237

uva7b$keep_remove

$keep
[1] "Defecate" "Itch"    

$remove
[1] "Urinate" "Tickle" 

Sample 7c

uva7c <- EGAnet::UVA(data = data7c, cut.off = 0.3)
uva7c

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.457
 Sneeze  Cough 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.280
 Urinate Defecate 0.254

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Hungry   Thirsty 0.242
  Vomit    Sneeze 0.224
  Heart Breathing 0.218

uva7c$keep_remove

$keep
[1] "Cough" "Itch" 

$remove
[1] "Sneeze" "Tickle"

Sample 8a

uva8a <- EGAnet::UVA(data = data8a, cut.off = 0.3)
uva8a

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.237
   Heart Breathing 0.218
  Hungry   Thirsty 0.213

uva8a$keep_remove

NULL

Sample 8b

uva8b <- EGAnet::UVA(data = data8b, cut.off = 0.3)
uva8b

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.277
   Heart Breathing 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.206

uva8b$keep_remove

NULL

Sample 9

uva9 <- EGAnet::UVA(data = data9, cut.off = 0.3)
uva9

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
  Tickle     Itch 0.379
    Wind     Burp 0.373
 Urinate Defecate 0.341

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.285
 Sneeze     Cough 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j  wto
 Hungry Thirsty 0.24

uva9$keep_remove

$keep
[1] "Defecate" "Wind"     "Itch"    

$remove
[1] "Urinate" "Burp"    "Tickle" 

Sample 10a

uva10a <- EGAnet::UVA(data = data10a, cut.off = 0.3)
uva10a

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.297

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.209

uva10a$keep_remove

NULL

Sample 10b

uva10b <- EGAnet::UVA(data = data10b, cut.off = 0.3)
uva10b

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva10b$keep_remove

NULL

Sample 10c

uva10c <- EGAnet::UVA(data = data10c, cut.off = 0.3)

Warning: Some variables did not belong to a dimension: Cough 

Use caution: These variables have been removed from the TEFI calculation

uva10c

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva10c$keep_remove

NULL

Sample 11a

uva11a <- EGAnet::UVA(data = data11a, cut.off = 0.3)
uva11a

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.330
   Heart Breathing 0.314
    Wind      Burp 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i  node_j   wto
 Muscles    Pain 0.279
  Tickle    Itch 0.277
  Hungry Thirsty 0.255

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
  Vomit Sneeze 0.201

uva11a$keep_remove

$keep
[1] "Heart"    "Defecate" "Burp"    

$remove
[1] "Breathing" "Urinate"   "Wind"     

Sample 11b

uva11b <- EGAnet::UVA(data = data11b, cut.off = 0.3)
uva11b

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
 Hungry Breathing 0.252

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j   wto
   Wind    Burp 0.214
   Temp Muscles 0.202

uva11b$keep_remove

NULL

Summary

rez_uva <- c(uva1a$keep_remove$remove,
  uva1b$keep_remove$remove,
  uva2$keep_remove$remove,
  uva3$keep_remove$remove,
  uva4$keep_remove$remove,
  uva5$keep_remove$remove,
  uva6$keep_remove$remove,
  uva7a$keep_remove$remove,
  uva7b$keep_remove$remove,
  uva7c$keep_remove$remove,
  uva8a$keep_remove$remove,
  uva8b$keep_remove$remove,
  uva9$keep_remove$remove,
  uva10a$keep_remove$remove,
  uva10b$keep_remove$remove,
  uva10c$keep_remove$remove,
  uva11a$keep_remove$remove,
  uva11b$keep_remove$remove)

sort(table(rez_uva), decreasing = TRUE) / length(rez_uva) 

rez_uva
    Tickle      Cough    Urinate  Breathing       Burp     Sneeze       Wind 
0.30434783 0.13043478 0.13043478 0.08695652 0.08695652 0.08695652 0.08695652 
  Defecate       Itch 
0.04347826 0.04347826 

Structure Refinement

Discarded the following items: - Tickle: not present in the same dataset and consistently flagged as redundant in UVA analysis.

data_all <- select(data_all, -Tickle)
data1a <- select(data1a, -Tickle)
data1b <- select(data1b, -Tickle)
data2 <- select(data2, -Tickle)
data3 <- select(data3, -Tickle)
data4 <- select(data4, -Tickle)
data5 <- select(data5, -Tickle)
data6 <- select(data6, -Tickle)
data7a <- select(data7a, -Tickle)
data7b <- select(data7b, -Tickle)
data7c <- select(data7c, -Tickle)
# data8a <- select(data8a, -Tickle)
# data8b <- select(data8b, -Tickle)
data9 <- select(data9, -Tickle)
data10a <- select(data10a, -Tickle)
data10b <- select(data10b, -Tickle)
data10c <- select(data10c, -Tickle)
data11a <- select(data11a, -Tickle)
data11b <- select(data11b, -Tickle)

Initial

colors <- c("Heart" = "#F44336", "Breathing" = "#FF5722",
            "Hungry" = "#FF9800", "Thirsty" = "#FFC107",
            "Burp" = "#4CAF50", "Wind" = "#009688",
            "Urinate" = "#63824C", "Defecate" = "#795548",
            "Cough" = "#8BC34A", "Sneeze" = "#CDDC39", 
            "Bruise" = "#673AB7", "Blood Sugar" = "#3F51B5",
            "Muscles" = "#2196F3", "Pain" = "#00BCD4",
            "Sex arousal" = "#FF4081", "Temp" = "#9C27B0",
            "Vomit" = "#76FF03", "Taste" = "#00E676",
            "Affective touch" = "#FF1744", "Itch" = "#D500F9")

HCA

make_hclust <- function(data) {
  rez_pvclust <- pvclust::pvclust(data,
                 method.hclust = "complete",
                 method.dist = "correlation",
                 nboot = 1000, quiet=TRUE, parallel=TRUE)
  
  
  dendrogram <- as.dist(1 - cor(data, use = "pairwise.complete.obs")) |> 
  hclust(method = "complete") |> 
  tidygraph::as_tbl_graph() 


  # Process Nodes
  nodes <- as.list(dendrogram)$nodes |> 
    mutate(
      Size = ifelse(label != "", 10, NA),
      Item = str_replace(label, "_", " "),
      idx = 1:nrow(as.list(dendrogram)$nodes))
  
  # Central node 
  nodes[nodes$height == max(nodes$height), c("Item", "Size")] <- data.frame(Item="Central", Size=15)
  
  # Process Edges
  edges <- as.list(dendrogram)$edges
  edges$linewidth = datawizard::rescale(nodes[edges$from, ]$height, to = c(0.1, 1))
  
  p <- tidygraph::tbl_graph(nodes = nodes, edges = edges) |> 
    ggraph(layout = "dendrogram", circular = TRUE) +
    # geom_edge_diagonal(strength = 0.7, linewidth = 1) +
    geom_edge_elbow2(aes(edge_width=linewidth), color="#212121") +
    geom_node_point(aes(filter=Item %in% c("Central"), size = Size), color="#212121") +
    # geom_node_point(aes(filter=Group %in% c("Visceroception", "Awareness", "Deficit"), color=Group, size = Size)) +
    geom_node_text(aes(
      label = ifelse(Item != "", Item, NA),
      x = x * 1.10,
      y = y * 1.10,
      filter = label != "",
      angle = ifelse(
        x >= 0,
        asin(y) * 360 / 2 / pi,
        360 - asin(y) * 360 / 2 / pi
      ),
      hjust = ifelse(
        x >= 0, 0, 1
      ))) +
    geom_node_point(aes(filter = label != "", color=Item, size=Size), alpha = 1) +
    # geom_node_text(aes(label=idx)) +  # Debug
    scale_edge_width_continuous(range=c(1, 3), guide = "none") +
    scale_size_continuous(range=c(10, 10), guide = "none") +
    scale_color_manual(values = colors,
                       breaks = names(colors)) +
    ggtitle("Hierarchical Clustering Analysis (HCA)", subtitle = "Method = Correlation") +
    coord_equal(clip = "off", xlim = c(-1.25, 1.25), ylim = c(-1.25, 1.25)) +
    theme_void() + 
    # guides(color = guide_legend(override.aes = list(size = c(7.5, 5, 5, 3.5, 7.5, 5, 5, 3.5, 3.5, 3.5, 7.5, 5, 5, 3.5)))) +
    theme(legend.text = ggtext::element_markdown(),
          legend.title = element_blank(),
          legend.position = "none",
          # plot.title = element_blank(), 
          # plot.subtitle = element_blank()
          plot.title = element_text(face="bold"))

  
  list(pvclust = rez_pvclust, p = p)
}


rez_hclust <- make_hclust(data_all)
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

rez_hclust$p

EGA

make_ega <- function(data, cols = c("C1" = "#F44336", "C2" = "#FF9800", "C3" = "#795548", "C4" = "#CDDC39", "C5" = "#2196F3", "C6" = "#009688", "C7" = "#9C27B0")) {
  rez <- EGAnet::bootEGA(
    data = data,
    seed = 123,
    model = "glasso",
    algorithm = "leiden",
    EGA.type = "hierEGA",
    type = "resampling",
    plot.itemStability = FALSE,
    verbose = FALSE, 
    allow.singleton = TRUE
  )
  
  table <- EGAnet::net.loads(rez$EGA$lower_order)$std |> 
    as.data.frame() |> 
    rownames_to_column("Item") |>
    gt::gt() |>
    gt::tab_header(title = "EGA Loadings") |>
    gt::data_color(
      columns = -Item,
      method = "numeric",
      colors = scales::col_numeric(
        palette = c("red", "white", "green"),
        domain = c(-1, 0, 1)
      )) |> 
    gt::fmt_auto()
  
    
  nodes <- rez$EGA$lower_order$dim.variables |> 
    rename(name = items) |> 
    mutate(dimension = paste0("C", dimension),
           size = apply(EGAnet::net.loads(rez$EGA$lower_order)$std, 1, max)[name])
  
  loadings <- rez$EGA$lower_order$network |> 
    as.data.frame() |> 
    rownames_to_column("to") |>
    pivot_longer(-c(to), names_to = "from", values_to = "Loading") |> 
    filter(Loading > quantile(Loading, 1/3)) 
  
  g <- tidygraph::tbl_graph(nodes = nodes, edges = loadings, directed = FALSE) |> 
    mutate(name = str_replace(name, "_", " "),
           name = ifelse(name == "Sex arousal", "Sex", name)) 
  
  set.seed(1)
  layout <- ggraph::create_layout(g, layout = "fr", weights = abs(loadings$Loading))
  
  xrange <- max(layout$x) - min(layout$x)
  yrange <- max(layout$y) - min(layout$y)
  xmin <- min(layout$x) - xrange * 0.05
  xmax <- max(layout$x) + xrange * 0.05
  ymin <- min(layout$y) - yrange * 0.05
  ymax <- max(layout$y) + yrange * 0.05
  
  p <- layout |> 
    ggraph() +
    geom_edge_bend(aes(edge_width = Loading, edge_alpha = Loading), color = "#212121",
                   strength = 0.3) +
    geom_node_point(aes(size = size, color = dimension), alpha = 0.95) +
    geom_node_text(aes(label = name), size = 3, color = "white", fontface = "bold") +
    scale_size_continuous(range = c(23, 28)) +
    scale_edge_width(range = c(0.3, 4)) +
    scale_edge_alpha(range = c(0.1, 0.9)) +
    scale_color_manual(values = cols) +
    labs(title = "Exploratory Graph Analysis (EGA)", subtitle = "Method = Leiden") +
    theme_void() +
    theme(legend.position = "none",
          plot.title = element_text(face="bold")) +
    coord_cartesian(xlim =c(xmin, xmax), ylim = c(ymin, ymax)) 
    
  list(rez = rez, table = table, p = p)
}

rez_ega <- make_ega(data_all, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]]))
stab <- EGAnet::itemStability(rez_ega$rez)

rez_ega$table

EGA Loadings
Item	1	2	3	4	5	6	7
Heart	0.524	0.038	0.009	0.067	0.038	0.019	0.114
Breathing	0.524	0.193	0.046	0.185	0.055	0.002	0
Hungry	0.133	0.522	0.07	0.124	0	0.003	0.037
Thirsty	0.095	0.522	0.124	0.105	0.009	0.011	0.01
Urinate	0.021	0.179	0.563	0.15	0.038	0.015	0
Defecate	0.046	0.057	0.563	0.156	0.113	0.042	0
Pain	0.037	0.109	0.055	0.387	0.046	0	0.138
Temp	0.026	0.089	0.086	0.365	0.127	0.028	0.007
Muscles	0.103	0.024	0.039	0.359	0.058	0.089	0.127
Sex_arousal	0.073	0.065	0.071	0.298	0.079	0.078	0
Affective_touch	0.037	0	0.034	0.276	0.05	0	0.207
Taste	0.07	0.03	0.062	0.225	0.116	0.039	0.055
Sneeze	0.03	0	0.051	0.108	0.569	0.081	0.047
Cough	0.038	0.012	0.012	0.136	0.446	0.162	0.086
Vomit	0.047	0	0.094	0.188	0.293	0.043	0.042
Wind	0.008	0.009	0.055	0.116	0.104	0.572	0.059
Burp	0.019	0.01	0.005	0.099	0.184	0.572	0.111
Bruise	0.03	0	0	0.186	0.048	0.081	0.44
Blood_Sugar	0.067	0.049	0	0.077	0.02	0.025	0.387
Itch	0.019	0	0	0.138	0.076	0.033	0.298

rez_ega$p

EFA

make_efa <- function(data, n = 3, cols = c("F1" = "#009688", "F2" = "#795548", "F3" = "red")) {
  n_fac <- n_factors(data)

  rez_efa <- factor_analysis(data, n = n, sort = TRUE, rotation = "oblimin")
  
  # Exctract loadings
  loadings <- attributes(rez_efa)$loadings_long |> 
    select(from = Component, to=Variable, Loading) |> 
    mutate(Type = "Loading",
           to = str_replace(to, "_", " "),
           from = str_replace(from, "MR", "F")) 
  
  # Reorder factors
  fct_order <- c("Pain", "Muscles", "Blood Sugar", "Bruise",
                 "Affective touch", "Sex arousal", "Temp", "Itch", "Tickle", 
                 "Thirsty", "Hungry", "Breathing", "Heart", "Defecate", "Urinate",
                 "Sneeze", "Cough",  "Wind", "Burp", "Vomit", "Taste")
  loadings <- loadings |> 
    mutate(to = fct_relevel(to, fct_order[fct_order %in% unique(loadings$to)])) |> 
    arrange(to) |> 
    mutate(to = as.character(to))

  # Get correlations between factors
  phi <- attributes(rez_efa)$model$Phi
  phi[upper.tri(phi, diag = TRUE)] <- NA
  phi <- phi |> 
    as.data.frame() |> 
    rownames_to_column("from") %>%
    pivot_longer(-from, names_to = "to", values_to = "Loading") |> 
    filter(!is.na(Loading)) |>
    mutate(Type = "Correlation",
           from = str_replace(from, "MR", "F"),
           to = str_replace(to, "MR", "F")) 
  
  # Reverse so that arcs are on the same side
  phi[paste(phi$from, phi$to) %in% c("F3 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F3")
  phi[paste(phi$from, phi$to) %in% c("F2 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F2")
  phi[paste(phi$from, phi$to) %in% c("F3 F2"), c("from", "to")] <- data.frame(from = "F2", to = "F3")
  phi[paste(phi$from, phi$to) %in% c("F4 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F4")
  phi[paste(phi$from, phi$to) %in% c("F4 F2"), c("from", "to")] <- data.frame(from = "F2", to = "F4")
  phi[paste(phi$from, phi$to) %in% c("F4 F3"), c("from", "to")] <- data.frame(from = "F3", to = "F4")
  
  loadings <- rbind(loadings, phi)
  
  # Make layout
  factor_names <- unique(loadings$from[loadings$Type == "Loading"])
  item_names <- unique(loadings$to[loadings$Type == "Loading"])
  
  nodes <- tibble(name = unique(c(loadings$from, loadings$to))) |>
    mutate(type = case_when(
      name %in% factor_names ~ "Factor",
      name %in% item_names ~ "Item",
      TRUE ~ "Other"
    ))
  
  y_position <- seq(0.8, 0.2, length.out = length(factor_names))
  
  layout <- nodes |>
    mutate(x = ifelse(type == "Item", 0, 1),
           y = case_when(
            name == "F1" ~ y_position[1],
            name == "F2" ~ y_position[2],
            name == "F3" ~ y_position[3],
            name == "F4" ~ y_position[4],
            name == "F5" ~ y_position[5],
            name == "F6" ~ y_position[6],
            .default = NA
          ))
  layout[layout$type == "Item", "y"] <- seq(0, 1, length.out = sum(layout$type == "Item"))
  
  
  g <- tidygraph::as_tbl_graph(loadings, directed = TRUE) |>
    left_join(layout, by = c("name"))
  
  p_efa <- ggraph(g, layout = as.data.frame(layout)[c("x", "y")]) + 
    geom_edge_link(
      aes(
        edge_width = abs(Loading),
        edge_alpha = abs(Loading),
        # edge_color = Loading,
        filter = Type == "Loading"
        # label = insight::format_value(Loading, lead_zero = FALSE, zap_small = TRUE)
      ),
      color = "#212121",
      arrow = arrow(length = unit(2.5, "mm"), type = "closed"),
      end_cap = circle(5, 'mm')
    ) +
      geom_edge_arc(
      aes(
        edge_width = abs(Loading),
        filter = Type == "Correlation",
        label = insight::format_value(Loading, lead_zero = FALSE, zap_small = TRUE)
      ),
      edge_color = "forestgreen",
      angle_calc = "along",
      label_dodge = unit(2.5, "mm"),
      force_flip = TRUE,
      strength = 0.2  # controls curvature
    ) +
    geom_node_label(data=filter(layout, type == "Item"), 
                    aes(fill = name), 
                    label = paste(rep(" ", 20), collapse = ""),
                    size = 3, 
                    hjust = 0.5, 
                    nudge_x = -0.05,
                    label.r = unit(0.2, "lines"),
                    label.padding = unit(0.5, "lines"),
                    label.size = 0) +
    geom_node_label(data=filter(layout, type == "Item"), 
                    aes(label = name), 
                    fill = NA,
                    size = 3, 
                    hjust = 0.5, 
                    vjust = 0.5,
                    color = "white",
                    fontface = "bold",
                    nudge_x = -0.05,
                    label.size = 0) +
    geom_node_point(data=filter(layout, type == "Factor"), 
                    aes(color = name), 
                    size = 12) +
    geom_node_text(data=filter(layout, type == "Factor"), 
                   aes(label = name), 
                   size = 4, 
                   color = "white", 
                   fontface = "bold") +
    scale_edge_width(range = c(0.3, 2)) +
    scale_edge_alpha(range = c(0.1, 1)) +
    scale_fill_manual(values = colors) +
    scale_color_manual(values = cols) +
    labs(title = "Exploratory Factor Analysis (EFA)", subtitle = "Method = Oblimin") +
    theme_void() +
    theme(legend.position = "none",
          plot.title = element_text(face="bold")) +
    coord_cartesian(xlim = c(-0.08, 1.08)) 

  list(n_fac=n_fac, efa=rez_efa, p=p_efa)
}


rez_efa <- make_efa(data_all, n = 3, cols = c("F1" = colors[["Wind"]], 
                                              "F2" = colors[["Defecate"]], 
                                              "F3" = colors[["Heart"]]))

Loading required namespace: GPArotation

plot(rez_efa$n_fac)

rez_efa$p

Final

We discarded:

• Taste: Lone item + unstable
• Affective_touch: Cross-loadings + unstable
• Vomit: Less strongly associated
• Itch: Less strongly associated, did not form consistent cluster
• More controversial: Temp and Sex_arousal: similar cluster yet less reliable

data_allf <- select(data_all, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data1af <- select(data1a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data1bf <- select(data1b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data2f <- select(data2, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data3f <- select(data3, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data4f <- select(data4, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data5f <- select(data5, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data6f <- select(data6, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data7af <- select(data7a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data7bf <- select(data7b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data7cf <- select(data7c, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data8af <- select(data8a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data8bf <- select(data8b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data9f <- select(data9, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data10af <- select(data10a, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data10bf <- select(data10b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data10cf <- select(data10c, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data11af <- select(data11a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data11bf <- select(data11b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)

HCA

rez_hclust <- make_hclust(data_allf)

plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

rez_hclust$p

EGA

rez_ega <- make_ega(data_allf, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))

stab <- EGAnet::itemStability(rez_ega$rez)

rez_ega$table

EGA Loadings
Item	1	2	3	4	5	6	7
Heart	0.528	0.042	0.017	0.041	0.026	0.042	0.134
Breathing	0.528	0.203	0.065	0.054	0.01	0.17	0.01
Hungry	0.144	0.525	0.083	0.008	0.013	0.096	0.053
Thirsty	0.098	0.525	0.128	0.018	0.018	0.104	0.017
Urinate	0.032	0.192	0.572	0.043	0.025	0.116	0
Defecate	0.069	0.072	0.572	0.087	0.062	0.115	0
Sneeze	0.057	0.009	0.089	0.563	0.104	0.087	0.03
Cough	0.054	0.021	0.036	0.563	0.176	0.12	0.078
Wind	0.014	0.024	0.071	0.117	0.576	0.074	0.067
Burp	0.032	0.015	0.016	0.179	0.576	0.098	0.127
Muscles	0.125	0.035	0.067	0.082	0.107	0.494	0.16
Pain	0.054	0.135	0.09	0.066	0.008	0.494	0.198
Bruise	0.042	0	0	0.051	0.092	0.281	0.488
Blood_Sugar	0.075	0.057	0	0.023	0.035	0.065	0.488

rez_ega$p

EFA

rez_efa <- make_efa(data_allf, n = 3, cols = 
                      c("F1" = colors[["Wind"]], "F2" = colors[["Defecate"]], 
                        "F3" = colors[["Heart"]]))
plot(rez_efa$n_fac)

rez_efa$p

CFA

make_cfa <- function(data) {
  
  # Main models
  m_g <- '
IAS =~ Burp + Wind + Cough + Sneeze + Hungry + Thirsty  + Urinate + Defecate + Heart + Breathing + Pain + Muscles 
  '
  m_ega <- '
    CoughSneeze =~ Cough + Sneeze
    UrinateDefecate =~ Urinate + Defecate
    BruiseBloodsugar =~ Bruise + Blood_Sugar
    WindBurp =~ Wind + Burp
    HeartBreath =~ Heart + Breathing
    MusclesPain =~ Muscles + Pain
    HungryThirsty =~ Hungry + Thirsty
  '
  
  m_efa <- '
    UrinateDefecate =~ Urinate + Defecate
    Sudden =~ Cough + Sneeze + Wind + Burp
    Homeostasis =~ Heart + Breathing + Muscles + Pain + Hungry + Thirsty + Bruise + Blood_Sugar
  '
    
    
  m_hclust <- '
    HeartBreath =~ Heart + Breathing
    BruiseBloodsugar =~ Bruise + Blood_Sugar
    HungryThirsty =~ Hungry + Thirsty
    Sudden =~ Cough + Sneeze + Wind + Burp
    Unpleasant =~ Urinate + Defecate + Muscles + Pain
  '
  
  # Test bifactor models
  m_egag <- paste0(m_ega, "\nIAS =~ CoughSneeze + WindBurp + BruiseBloodsugar + UrinateDefecate + HeartBreath + MusclesPain")
  m_efag <- paste0(m_efa, "\nIAS =~ UrinateDefecate + Sudden + Homeostasis")
  m_hclustg <- paste0(m_hclust, "\nIAS =~ HeartBreath + BruiseBloodsugar + HungryThirsty + Sudden + Unpleasant")

  fit_g <- lavaan::cfa(m_g, data = data, std.lv = TRUE)
  fit_ega <- lavaan::cfa(m_ega, data = data, std.lv = TRUE)
  fit_egag <- lavaan::cfa(m_egag, data = data, std.lv = TRUE)
  fit_efa <- lavaan::cfa(m_efa, data = data, std.lv = TRUE)
  fit_efag <- lavaan::cfa(m_efag, data = data, std.lv = TRUE)
  fit_hclust <- lavaan::cfa(m_hclust, data = data, std.lv = TRUE)
  fit_hclustg <- lavaan::cfa(m_hclustg, data = data, std.lv = TRUE)
  
  # Test
  bf <- test_bf(fit_ega, fit_g, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg) |> 
    as.data.frame() |> 
    select(Name = Model, BF)
  
  # test <- lavaan::lavTestLRT(fit_g, fit_ega, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg)
  test <- compare_performance(fit_g, fit_ega, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg,
                      metrics = c("BIC", "RMSEA", "Chi2", "CFI")) |> 
    select(-BIC_wt, -Model) |> 
    merge(bf, by = "Name") |>
    arrange(BIC) |> 
    gt::gt() |> 
    gt::tab_header(title = "CFA Model Comparison") |>
    gt::data_color(
      columns = c(BIC, Chi2, RMSEA),
      method = "numeric",
      palette = c("green", "white")) |>
    gt::data_color(
      columns = c(CFI),
      method = "numeric",
      palette = c("white", "green")) |>
    gt::data_color(
      columns = c(BF),
      method = "numeric",
      palette = c("red", "yellow", "white", "yellow", "green"),
      domain = c(0, 0.3, 1, 3, 10000000)) |>
    gt::fmt_number(columns = c(RMSEA, Chi2, CFI), decimals = 3) 

  
  param  <- parameters::parameters(fit_ega)
  perf <- performance::performance(fit_ega)
  
  # Graph
  nodes <- tidySEM::get_nodes(fit_ega) |> 
    mutate(name = str_replace(name, "_", " "))
  edges <- tidySEM::get_edges(fit_ega) |> 
    mutate(from = str_replace(from, "_", " "),
           to = str_replace(to, "_", " "), 
           est_std = as.numeric(est_std)) |> 
    filter(from != to)
  
  g <- tidygraph::tbl_graph(nodes = nodes, edges = edges, directed = TRUE) 
  set.seed(1)
  layout <- ggraph::create_layout(g, layout = "fr", weights = est_std)
  
  p_cfa <- ggraph(layout) +
    geom_edge_bend(aes(edge_width = abs(est_std), edge_alpha = abs(est_std),
                       label = format_value(est_std, zap_small = TRUE, lead_zero = FALSE),
                       filter = arrow == "none"), 
                   color = "forestgreen",
                   angle_calc = "along", 
                   check_overlap = TRUE,
                   label_dodge = unit(0, "mm"),
                   label_size = unit(3, "mm"),
                   strength = 0.3) +
    geom_edge_link(aes(edge_width = abs(est_std), edge_color = abs(est_std), 
                       label = format_value(est_std, zap_small = TRUE, lead_zero = FALSE),
                      filter = arrow != "none"), 
                   angle_calc = "along", 
                   label_dodge = unit(3, "mm"),
                   arrow = arrow(length = unit(2.5, "mm"), type = "open", ends = "first"), 
                   lineend = "round",
                   linejoin = "bevel",
                   start_cap = circle(9, 'mm')) +
    geom_node_point(aes(color = name, shape=shape), alpha = 0.97, size = 20) +
    geom_node_text(data = filter(layout, shape == "rect"), aes(label = name), 
                   size = 3, color = "white", fontface = "bold") +
    scale_size_continuous(range = c(10, 20)) +
    scale_edge_width(range = c(0.3, 4)) +
    scale_edge_alpha(range = c(0.1, 0.9)) +
    scale_edge_color_gradient(low = "#EEEEEE", high = "#212121") +
    scale_color_manual(values = c(colors, 
                               "CoughSneeze" = "#CDDC39", 
                               "UrinateDefecate" = "#795548", "BruiseBloodsugar" = "#673AB7",
                               "WindBurp" = "#009688", "HeartBreath" = "#F44336",
                               "MusclesPain" = "#2196F3", "HungryThirsty" = "#FF9800")) +
    scale_shape_manual(values = c("oval" = 16, "rect" = 15)) +
    labs(title = "Confirmatory Factor Analysis (CFA)", subtitle = "Method = Maximum Likelihood") +
    theme_void() +
    theme(legend.position = "none",
          plot.title = element_text(face="bold"))
  p_cfa
  
  list(test = test, param = param, perf = perf, p=p_cfa)
}


rez_cfa0 <- make_cfa(data_allf)

Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.

Warning: When comparing models, please note that probably not all models were fit
  from same data.

When comparing models, please note that probably not all models were fit
  from same data.

rez_cfa0$test

CFA Model Comparison
Name	BIC	RMSEA	Chi2	CFI	BF
fit_ega	-158722.7	0.036	2,290.274	0.984	NA
fit_egag	-154419.8	0.055	6,738.245	0.952	0
fit_hclust	-148000.0	0.079	13,126.977	0.905	0
fit_hclustg	-146829.7	0.079	14,349.064	0.896	0
fit_g	-145655.5	0.125	26,541.561	0.783	0
fit_efa	-144933.2	0.083	16,266.333	0.883	0
fit_efag	-144933.2	0.083	16,266.333	0.883	0

rez_cfa0$perf

# Indices of model performance

Chi2(56) | p (Chi2) | Baseline(91) | p (Baseline) |   GFI |  AGFI |   NFI
-------------------------------------------------------------------------
2290.274 |   < .001 |    1.380e+05 |       < .001 | 0.990 | 0.981 | 0.983

Chi2(56) |  NNFI |   CFI | RMSEA |      RMSEA  CI | p (RMSEA) |   RMR |  SRMR
-----------------------------------------------------------------------------
2290.274 | 0.974 | 0.984 | 0.036 | [0.034, 0.037] |    > .999 | 0.001 | 0.020

Chi2(56) |   RFI |  PNFI |   IFI |   RNI | Loglikelihood |        AIC
---------------------------------------------------------------------
2290.274 | 0.973 | 0.605 | 0.984 | 0.984 |     79615.220 | -1.591e+05

Chi2(56) |        BIC | BIC_adjusted
------------------------------------
2290.274 | -1.587e+05 |   -1.589e+05

rez_cfa0$param

# Loading

Link                            | Coefficient |       SE |       95% CI
-----------------------------------------------------------------------
CoughSneeze =~ Cough            |        0.17 | 1.16e-03 | [0.17, 0.17]
CoughSneeze =~ Sneeze           |        0.15 | 1.10e-03 | [0.15, 0.15]
UrinateDefecate =~ Urinate      |        0.17 | 1.23e-03 | [0.17, 0.17]
UrinateDefecate =~ Defecate     |        0.16 | 1.19e-03 | [0.16, 0.16]
BruiseBloodsugar =~ Bruise      |        0.22 | 2.11e-03 | [0.21, 0.22]
BruiseBloodsugar =~ Blood_Sugar |        0.16 | 1.94e-03 | [0.15, 0.16]
WindBurp =~ Wind                |        0.18 | 1.29e-03 | [0.18, 0.18]
WindBurp =~ Burp                |        0.20 | 1.32e-03 | [0.20, 0.20]
HeartBreath =~ Heart            |        0.14 | 1.49e-03 | [0.14, 0.14]
HeartBreath =~ Breathing        |        0.17 | 1.44e-03 | [0.17, 0.18]
MusclesPain =~ Muscles          |        0.15 | 1.19e-03 | [0.15, 0.15]
MusclesPain =~ Pain             |        0.15 | 1.23e-03 | [0.14, 0.15]
HungryThirsty =~ Hungry         |        0.18 | 1.62e-03 | [0.18, 0.18]
HungryThirsty =~ Thirsty        |        0.18 | 1.58e-03 | [0.18, 0.19]

Link                            |      z |      p
-------------------------------------------------
CoughSneeze =~ Cough            | 147.35 | < .001
CoughSneeze =~ Sneeze           | 136.55 | < .001
UrinateDefecate =~ Urinate      | 135.80 | < .001
UrinateDefecate =~ Defecate     | 136.25 | < .001
BruiseBloodsugar =~ Bruise      | 102.29 | < .001
BruiseBloodsugar =~ Blood_Sugar |  80.90 | < .001
WindBurp =~ Wind                | 141.15 | < .001
WindBurp =~ Burp                | 150.94 | < .001
HeartBreath =~ Heart            |  95.60 | < .001
HeartBreath =~ Breathing        | 119.51 | < .001
MusclesPain =~ Muscles          | 124.42 | < .001
MusclesPain =~ Pain             | 119.27 | < .001
HungryThirsty =~ Hungry         | 112.05 | < .001
HungryThirsty =~ Thirsty        | 115.23 | < .001

# Correlation

Link                                | Coefficient |       SE |       95% CI
---------------------------------------------------------------------------
CoughSneeze ~~ UrinateDefecate      |        0.56 | 5.94e-03 | [0.55, 0.57]
CoughSneeze ~~ BruiseBloodsugar     |        0.55 | 7.32e-03 | [0.53, 0.56]
CoughSneeze ~~ WindBurp             |        0.72 | 4.80e-03 | [0.71, 0.73]
CoughSneeze ~~ HeartBreath          |        0.52 | 6.73e-03 | [0.51, 0.53]
CoughSneeze ~~ MusclesPain          |        0.67 | 5.99e-03 | [0.65, 0.68]
CoughSneeze ~~ HungryThirsty        |        0.46 | 7.02e-03 | [0.45, 0.47]
UrinateDefecate ~~ BruiseBloodsugar |        0.40 | 7.83e-03 | [0.38, 0.41]
UrinateDefecate ~~ WindBurp         |        0.51 | 6.10e-03 | [0.50, 0.52]
UrinateDefecate ~~ HeartBreath      |        0.54 | 6.71e-03 | [0.52, 0.55]
UrinateDefecate ~~ MusclesPain      |        0.65 | 6.14e-03 | [0.64, 0.66]
UrinateDefecate ~~ HungryThirsty    |        0.65 | 6.15e-03 | [0.64, 0.66]
BruiseBloodsugar ~~ WindBurp        |        0.58 | 7.09e-03 | [0.57, 0.59]
BruiseBloodsugar ~~ HeartBreath     |        0.50 | 8.08e-03 | [0.48, 0.51]
BruiseBloodsugar ~~ MusclesPain     |        0.74 | 7.37e-03 | [0.73, 0.76]
BruiseBloodsugar ~~ HungryThirsty   |        0.45 | 8.30e-03 | [0.43, 0.47]
WindBurp ~~ HeartBreath             |        0.46 | 6.86e-03 | [0.44, 0.47]
WindBurp ~~ MusclesPain             |        0.63 | 6.03e-03 | [0.62, 0.65]
WindBurp ~~ HungryThirsty           |        0.43 | 7.01e-03 | [0.42, 0.45]
HeartBreath ~~ MusclesPain          |        0.64 | 6.88e-03 | [0.63, 0.66]
HeartBreath ~~ HungryThirsty        |        0.64 | 6.90e-03 | [0.63, 0.66]
MusclesPain ~~ HungryThirsty        |        0.63 | 6.90e-03 | [0.62, 0.65]

Link                                |      z |      p
-----------------------------------------------------
CoughSneeze ~~ UrinateDefecate      |  94.23 | < .001
CoughSneeze ~~ BruiseBloodsugar     |  74.98 | < .001
CoughSneeze ~~ WindBurp             | 150.87 | < .001
CoughSneeze ~~ HeartBreath          |  77.12 | < .001
CoughSneeze ~~ MusclesPain          | 111.01 | < .001
CoughSneeze ~~ HungryThirsty        |  65.50 | < .001
UrinateDefecate ~~ BruiseBloodsugar |  51.02 | < .001
UrinateDefecate ~~ WindBurp         |  83.50 | < .001
UrinateDefecate ~~ HeartBreath      |  79.79 | < .001
UrinateDefecate ~~ MusclesPain      | 105.72 | < .001
UrinateDefecate ~~ HungryThirsty    | 105.96 | < .001
BruiseBloodsugar ~~ WindBurp        |  81.97 | < .001
BruiseBloodsugar ~~ HeartBreath     |  61.56 | < .001
BruiseBloodsugar ~~ MusclesPain     | 100.55 | < .001
BruiseBloodsugar ~~ HungryThirsty   |  54.31 | < .001
WindBurp ~~ HeartBreath             |  66.62 | < .001
WindBurp ~~ MusclesPain             | 105.09 | < .001
WindBurp ~~ HungryThirsty           |  61.95 | < .001
HeartBreath ~~ MusclesPain          |  93.45 | < .001
HeartBreath ~~ HungryThirsty        |  93.47 | < .001
MusclesPain ~~ HungryThirsty        |  91.62 | < .001

rez_cfa0$p

Figure

The structural analysis seem to converge on the idea of small clusters (of pairs or triplets) that are potentially inter-related parts of larger clusters (although with unstable associations). Best way to assess the values of this new granular structure is to assess its predictive value against convergent measures, and see if it’s superior to a unique score (-> Study 2).

fig2 <- (rez_hclust$p | rez_ega$p) / (rez_efa$p | rez_cfa0$p) 
ggsave("figures/Figure2.png", fig2, width=14, height=14, dpi=300, bg="white")

Structure Invariance

Sample 1a

Initial

Summary

rez_hclust <- make_hclust(data1a)
rez_ega <- make_ega(data1a, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data1a, n = 4, cols = c("F1" = colors[["Wind"]], "F2" = colors[["Defecate"]], 
                                            "F3" = colors[["Heart"]], "F4" = colors[["Itch"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Details

plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

stab <- EGAnet::itemStability(rez_ega$rez)

rez_ega$table

EGA Loadings
Item	1	2	3	4	5	6	7	8
Heart	0.504	0.125	0.046	0	0.093	0.114	0	0
Breathing	0.504	0.14	0.032	0.098	0.158	0.123	0	0.036
Hungry	0.2	0.511	0.12	0.004	0.028	0.064	0	0.034
Thirsty	0.074	0.511	0.166	0	0.068	0.13	0.005	0.053
Defecate	0.044	0.093	0.632	0.107	0.077	0.036	0.059	−0.045
Urinate	0.026	0.253	0.421	0.035	0.167	0.035	0.046	−0.02
Taste	0.032	0.018	0.21	0.105	0.18	0.149	0.082	0.028
Sneeze	0	0.006	0.078	0.676	0.094	0.001	0.092	0.044
Cough	0.053	0	0.014	0.377	0.11	0	0.171	0.121
Vomit	0.088	0	0.177	0.299	0.066	0.067	0	0.054
Muscles	0.084	0.022	0.13	0.03	0.454	0.032	0.029	0.093
Pain	0.126	0.032	0.133	0.067	0.438	0.195	0	0.143
Temp	0.086	0.056	0.117	0.125	0.264	0.245	0	0.038
Sex_arousal	0.103	0.042	0.055	0.006	0.069	0.435	0.119	0.007
Affective_touch	0.059	0.086	0.058	0.027	0.205	0.435	0.011	0.152
Wind	0	0.007	0.189	0.051	0	0.261	0.571	0.058
Burp	0	0	0.007	0.203	0.034	0	0.571	0.105
Blood_Sugar	0.013	0	−0.02	0.04	0.033	0.033	0.015	0.455
Bruise	0	0.108	−0.044	0.08	0.193	0.044	0.068	0.401
Itch	0.034	0	0.027	0.076	0.074	0.224	0.07	0.388

plot(rez_efa$n_fac)

Final

Summary

rez_hclust <- make_hclust(data1af)
rez_ega <- make_ega(data1af, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data1af, n = 4, cols = c("F1" = colors[["Wind"]], "F2" = colors[["Defecate"]], 
                                            "F3" = colors[["Heart"]], "F4" = colors[["Itch"]]))
rez_cfa1a  <- make_cfa(data1af)

rez_hclust$p / rez_ega$p / rez_efa$p

Details

plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

stab <- EGAnet::itemStability(rez_ega$rez)

rez_ega$table

EGA Loadings
Item	1	2	3	4	5	6	7
Heart	0.511	0.128	0.036	0.009	0.01	0.108	0
Breathing	0.511	0.16	0.048	0.068	0.015	0.151	0.045
Hungry	0.202	0.514	0.105	0.015	0	0.042	0.048
Thirsty	0.091	0.514	0.178	0	0.018	0.041	0.089
Urinate	0.029	0.253	0.565	0.038	0.05	0.185	−0.017
Defecate	0.08	0.108	0.565	0.098	0.087	0.09	−0.042
Sneeze	0.033	0.02	0.12	0.574	0.117	0.091	0
Cough	0.071	0	0.021	0.574	0.171	0.06	0.157
Wind	0.021	0.024	0.136	0.071	0.577	0.028	0.089
Burp	0.012	0	0.008	0.22	0.577	0.041	0.089
Muscles	0.091	0.042	0.118	0.058	0.04	0.523	0.125
Pain	0.183	0.044	0.107	0.061	0.014	0.523	0.18
Bruise	0.01	0.128	−0.031	0.059	0.093	0.222	0.502
Blood_Sugar	0.033	0	−0.013	0.052	0.032	0.053	0.502

plot(rez_efa$n_fac)

rez_cfa1a$test

CFA Model Comparison
Name	BIC	RMSEA	Chi2	CFI	BF
fit_ega	-1635.181	0.026	72.825	0.988	NA
fit_egag	-1613.162	0.059	180.403	0.924	1.654961e-05
fit_hclustg	-1594.961	0.065	210.828	0.905	1.846821e-09
fit_hclust	-1592.739	0.062	182.492	0.921	6.080789e-10
fit_efag	-1564.540	0.073	253.472	0.877	4.576430e-16
fit_efa	-1564.540	0.073	253.472	0.877	4.576430e-16
fit_g	-1538.631	0.106	325.704	0.784	1.082552e-21

Sample 1b

TODO.

Model Performance

rbind(
  mutate(rez_cfa0$perf, Sample = "All Samples"),
  mutate(rez_cfa1a$perf, Sample = "Sample 1a")
  # TODO
) |> 
  select(Sample, Chi2, RMSEA, CFI, "SRMR") 

# Indices of model performance

Sample      |     Chi2 | RMSEA |   CFI |  SRMR
----------------------------------------------
All Samples | 2290.274 | 0.036 | 0.984 | 0.020
Sample 1a   |   72.825 | 0.026 | 0.988 | 0.029

Scores

Note the usage of descriptive factor names relating directly to the items to avoid abstraction.

# Add empirical variables
add_empirical <- function(data, sample = "Sample 1a") {
  x <- data.frame(
    Original = rowMeans(data),
    HungryThirsty = (data$Hungry + data$Thirsty) / 2,
    MusclesPain = (data$Muscles + data$Pain) / 2,
    WindBurp = (data$Wind + data$Burp) / 2,
    UrinateDefecate = (data$Urinate + data$Defecate) / 2,
    BreathingHeart = (data$Heart + data$Breathing) / 2)
  
  if("Blood_Sugar" %in% names(data)) {
    x$BruiseBlood <- (data$Blood_Sugar + data$Bruise) / 2
  } else {
    x$BruiseBlood <- data$Bruise
  }
  
  if("Cough" %in% names(data)) {
    x$CoughSneeze <- (data$Cough + data$Sneeze) / 2
  } else {
    x$CoughSneeze <- data$Sneeze
  }
  x$Sample <- sample
  x
}

scores <- list(
  sample1a = add_empirical(data1a, "Sample 1a"),
  sample1b = add_empirical(data1b, "Sample 1b"),
  sample2 = add_empirical(data2, "Sample 2"),
  sample3 = add_empirical(data3, "Sample 3"),
  sample4 = add_empirical(data4, "Sample 4"),
  sample5 = add_empirical(data5, "Sample 5"),
  sample6 = add_empirical(data6, "Sample 6"),
  sample7a = add_empirical(data7a, "Sample 7a"),
  sample7b = add_empirical(data7b, "Sample 7b"),
  sample7c = add_empirical(data7c, "Sample 7c"),
  sample8a = add_empirical(data8a, "Sample 8a"),
  sample8b = add_empirical(data8b, "Sample 8b"),
  sample9 = add_empirical(data9, "Sample 9"),
  sample10a = add_empirical(data10a, "Sample 10a"),
  sample10b = add_empirical(data10b, "Sample 10b"),
  sample10c = add_empirical(data10c, "Sample 10c"),
  sample11a = add_empirical(data11a, "Sample 11a"),
  sample11b = add_empirical(data11b, "Sample 11b")
)

save(scores, file = "../data/scores.RData")